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NOTATION 


K  =  Number  of  blades  of  the  propeller 

k  =  Reduced  frequency 
n 

K  =  Nondimensional  pressure  coefficient,   K  =  *— 


P  P    4pf(RPM/60)2R2 

£  =  Blade  section  length 

N  =  Harmonic  number 

NT  =  Nondimensional  time,  NT  =  t  X  RPM/60 

P  =  Pitch  in  general 

Ap  =  Pressure  jump  in  general 

R  =  Propeller  radius 

r,  =  Hub  radius 
h 

RPM  =  Propeller  revolutions  per  minute 

s,r;n  =  Curvilinear  coordinate  system  defined  in  Fig.  2. 

T  =  Maximum  thickness  of  propeller  blade  section 

t  =  Time  variable 

u  =  Axial  perturbation  velocity  in  general 

u  =  Tangentical  perturbation  velocity  in  general 

V  =  Relative  flov;  past  a  blade  section 
o  r 

V  =  Ship  speed 

x,r,0  =  Cylindrical  coordinate  system  defined  in  Fig.  2 

x,y,z  =  Cartesian  coordinate  system  defined  in  Fig.  2 

P  =  Advance  angle 

T  -  Circulation  in  general 

6,  =  A  coordinate  of  tip  of  k'th  blade 


(vii) 


A  =  Advance  coefficient 

6  =6    coordinate  of  leading  edge 

Li 

0  =9    coordinate  of  trailing  edge 

ft  =  Propeller  angular  velocity 

0  =  Velocity  potential  in  general 

p  =  Fluid  density 

T  =  Blade  section  thickness 

l,p,cp  =  Dummy  cylindrical  coordinates 

£,T),£  =  Dummy  cartesian  coordinates 

Note:  Other  symbols  and  abbreviations  are  explained  whenever  used, 
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Chapter  I 
INTRODUCTION 

The  propeller  induced  vibratory  forces  acting  on  a  ship  have 
been  for  a  long  time  an  important  subject  of  experimental  and  theo- 
retical research.   Although  some  important  conclusions  have  been 
reached  based  on  reliable  experimental  results  [17],  it  is  always 
hard,  if  not  impossible,  to  support  such  conclusions  with  theoretical 
calculations.   One  reason  for  this  defficiency  is  that  for  a  long  time 
no  analytical  method  was  available  for  evaluating  the  propeller  loading 
as  a  three-dimensional  problem. 

The  initial  numerical  results  obtained  for  the  propeller-induced 
pressure  field  were  based  on  a  simple  lifting-line  representation  [10] 
which  was  found  to  predict  pressures  far  below  those  obtained  from 
measurement.   This  theory  was  later  refined  by  including  the  effect  of 
the  blade  thickness  and  this  brought  the  numerical  results  closer  to 
the  experimental  data  [9,  10,  11]. 

The  effects  on  the  pressure  field  due  to  nonuniform  inflow  were 
studied  by  Tsakonas  e_t  a_l.  [8].   They  treated  the  problem  as  a  Dirichlet 
problem  of  potential  theory  and  related  the  pressure  field  near  the 
propeller  to  the  trust-  and  torque-producing  forces.   These  forces 
were  computed  using  the  gust  theory  for  airfoil^   and  assuming  that 
the  flow  around  each  blade  was  two-dimensional.   Two  correction  factors 
were  added  to  compensate  for  cascade  and  finite  aspect  ratio  effects. 
Even  though  the  calculations  made  using  this  method  cannot  be  con- 
sidered very  accurate,  they  can  give  an  idea  of  the  relative  importance 
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of  the  unsteady  effects  arising  from  the  nonunif ormity  of  the  ship's 
wake.   This  can  be  seen  in  Fig.  1,  taken  from  [8].   It  is  apparent 
that  very  close  to  the  plane  of  the  propeller  the  increase  in  pressure 
due  to  nonuniformity  is  at  most  20  percent  of  the  total  pressure  but 
beyond  x/D  =0.4  the  percentage  change  becomes  much  greater.   The 
reason  is  that  at  such  distances  all  the  pressures  are  small,  but  the 
slow  rate  of  decay  of  the  components  arising  from  nonuniformity  makes 
them  more  important.   It  must  be  noted,  however,  that  the  pressure  due 
to  thickness  is  not  included  in  Fig.  1  so  that  the  relative  effect 
close  to  the  propeller  arising  from  nonuniform  inflow  is  of  the  order 
of  10  percent  for  the  particular  wake  used  in  these  calculations. 

A  truly  three-dimensional  theory  was  recently  used  by  Kerwin 
[18]  to  compute  the  field  point  velocities  induced  by  a  marine  pro- 
peller operating  in  a  steady  flow.   He  makes  use  of  the  vortex  theory 
to  represent  the  blade  loading  and  derives  a  scheme  to  compute  the  per- 
turbation velocities  on  specified  field  points  located  on  the  blades. 

In  the  present  work  we  will  extend  the  program  developed  by 
Professor  Kerwin  to  include  field  points  located  at  any  point  in  the 
fluid  around  the  propeller  and  to  take  into  account  the  unsteady  effects 
due  to  a  nonuniform  inflow  such  as  is  found  in  a  ship's  wake.   In 
addition,  a  scheme  will  be  developed  to  compute  the  velocity  potential 
at  the  field  points.   The  pressure  field  is  then  computed  by  sub- 
stituting the  perturbation  velocities  and  time  variation  of  the 
velocity  potential  into  Bernoulli's  equation. 
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The  usual  assumptions  of  an  incompressible,  frictionless  and 
irrotational  flow  are  held.   It  is  also  assumed  that  all  the  disturb- 
ances are  small,  which  requires  that  the  blades  be  thin  and  propeller 
loading  small.   The  effects  of  thickness  form  and  outline  of  the 
blades  and  the  form  of  the  loading  distribution,  both  radially  and 
over  the  chord  are  taken  into  account  by  this  theory. 

When  this  work  was  already  under  way,  the  Davidson  Laboratory 
published  a  report  by  Jacobs  e_t  aT.  [19]  which  also  treates  the  un- 
steady, propeller  induced,  pressure  field  using  a  lifting  surface 
theory.   However,  their  treatment  is  based  on  the  solution  of  a  Dirichlet 
problem  for  the  pressure  field.   It  seemed  of  interest  to  compare  the 
results  obtained  with  the  two  different  approaches  but  unfortunately 
were  were  not  able  to  obtain  from  D.L.  the  loading  distribution  for 
the  propellers  used  in  their  report. 
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Chapter  II 
FORMULATION  OF  THE  PROBLEM 

It  is  desired  to  find  the  pressure  at  any  point  in  the  fluid 
around  a  ship's  propeller.   Assuming  that  the  flow  is  incompressible, 
frictionless  and  irrotational,  the  problem  becomes  one  of  potential 
theory  which  in  principle  can  be  solved  by  finding  the  velocity 
potential  and  velocities  induced  by  the  propeller  and  substituting 
into  the  Bernoulli's  equation  to  solve  for  the  pressure. 

1.   Geometry  of  the  Propeller 

The  following  geometrical  considerations  are  adopted  from  [1]. 

Two  orthogonal  coordinate  systems  --  a  Cartesian  (x,y,z)  and  a 
cylindrical  (x,r,8)  --  are  associated  with  the  propeller  as  represented 
in  Fig.  2.   Both  systems  are  fixed  with  the  propeller  and  the  first 
blade  is  made  to  coincide  with  the  y  axis.   All  length  dimensions  are 
nondimensionalized  with  respect  to  the  radius  of  the  propeller,  R. 
The  equations  relating  the  two  coordinate  systems  are: 

2  2   - 

y  =  r  cos0      r  =  (y  +  z  )2 

(1) 
z  =  r  sin0      Q   =  tan"1  (z/y) 

In  order  to  relate  corresponding  points  on  each  of  the  K  blades, 
we  define  5,  as  the  6    coordinate  of  the  point  at  the  tip  of  the  k'th 
blade.   For  symmetrically  arranged  blades  it  is 


S^11^       k=l,2,...,K  (2) 
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FIG.    2   -  COORDINATE  SYSTEMS 
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FIG.  3  -  PROJECTED  BLADE  OUTLINE 
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This  is  illustrated  in  Fig.  3  which  also  shows  the  projected  blade 
outline,  defined  by  the  angular  coordinates  0  (r)  and  0  (r)  of  the 
leading  and  trailing  edges,  respectively. 

The  geometry  of  the  flow  past  a  blade  is  schematically  re- 
presented in  Fig.  h.      It  is  seen  that  the  oncoming  flow  V  (r)  forms 
an  angle  (3  with  the  yz  plane,  and  we  define 

Mr)  =  r  tanP(r)  (3) 

as  the  "advance  coefficient".   The  streamlines  are  helices  with 
pitch 

P(r)  =  27T  A(r)  -  (10 

It  is  useful  to  define  an  orthogonal  curvilinear  coordinate 
system  (s,r,n)  associated  with  the  helical  surfaces  swept  out  by 
the  undisturbed  flow  past  the  radial  lines  {x  =  0,  0  =  5,  ).   This 
system,  shown  in  Fig.  2,  has  unit  vectors: 

u  =i  sinp  -   sin(0  +5  )cosP +k  cos(0  +5,  )cosP 

S  K.  K 

u  =  cos(0+5,)     +  k  sin(0  +5,)  (5) 

u  =  i  cos{3  +   sin(0  +6,  )  sinp  -  k  cos(0  +5,  )  sinf3 
n  k  k 

It  is  assumed  that  the  blades  lie  on  those  helical  surfaces. 
This  assumption  is  in  accordance  with  the  linear  theory  for  lightly 
loaded  propellers,  however,  it  does  not  impose  any  limitation  to  the 
possibility  of  including  nonlinear  refinements  by  suitably  modifying 
the  advance  angle  ['!]. 


NO  ORDER  SLIPS  INVOLVED 
IN  TECHNICAL  PROCESSING 
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FIG.  4  -  APROACH  FLOW  VELOCITY  DIAGRAM 
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2 .   Expression  for  the  Pressure 

The  Bernoulli's  equation  for  unsteady  flow  is 

^         2     2 

p  =  pf  ^  +  i  pf (va  -  v  > 


where 


p   =  hydrodynamic  pressure 

pf  =  fluid  density 

<J>  =  velocity  potential 

V.  =  speed  of  free  stream  or  speed  of  advance 

V  =  resultant  velocity  at  a  point 

t   =  time  variable. 


Changing  to  a  frame  of  reference  moving  with  the  propeller 
blades  and  after  linearization  the  above  equation  becomes  [11] 

P  =  Pf  ^  +  PfUaVA  +  pf°r  Ut  (6) 

where 

u  ,u   =  axial  and  tangential  disturbance  velocities 
a   t 

Q  =  angular  velocity  of  the  propeller. 

3 .   Induced  Velocities 

In  order  to  find  the  velocities  induced  by  the  propeller,  we 
make  use  of  the  vortex  representation  described,  for  instance,  in 
[1]  and  [2].   The  following  is  an  adaptation  to  the  unsteady  case 
of  the  relations  derived  in  [1]. 
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The  inflow  velocity  and  verticity  strength  are  expressed  as 
complex  Fourier  series  of  the  angular  position  of  the  rotating  co- 
ordinate system  with  respect  to  a  set  of  axes  fixed  in  the  ship, 

VQ(p,cpt)  =  Vo(p)-e 

iNfit 


r(p,cp,t)  =  r(p,cp)-e' 


where 


V  =  relative  flow  velocity  past  the  blade  section 
o 

Y  =   vorticity  strength,  in  general 
N   =  harmonic  number 

n  =  angular  velocity  of  the  propeller 

(i)  Velocity  Induced  by  Source  Distribution 

The  velocity  induced  at  a  point  (x,r,0)  by  sources  distributed 
over  K  blades  is 

K        iN(fic-5k)    O1    pT 

Ua      (x,r,9,t)    f,  kr  J     1      V<>>  " 

k=1  rh     8L 

.a-«p.  |i .  r»-MP)-9i-rp2+^(p)i*  ^  dp  (7) 

D 

iN(J2t-6k)     i       e 
K  P~    O    T 

u(s)(x,r,e,t)  -E ;^Tf J    J        vo(P) 

k=1  ^h  eL 

a  p    sin(cp+5k-e).[p2-r*2(p)r 

•e     Y-    v-    •    3 dcp    dp  (8) 

os  Do 
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where 

D  =  [(x  -A(p)'cp)2+  r2+  p2-  2rp  cos(cp+5k-  0)]2 

X  =  blade  thickness 
The  superscript  s  denotes  a  source  velocity. 

(ii)  Velocity  Induced  by  Vortex  Distribution 
We  separate  the  components  of  the  induced  velocities  due  to 
the  bound  vorticity  in  the  blade  (superscript  b)  ,   radial  vorticity 
in  the  wave  (sup.  r),  streamwise  vorticity  in  the  blade  (sup.  f) 
and  streamwise  vorticity  in  the  wake  (sup.  w) . 

K    iN(ftt-5k)  p1  p0T 
u<b)(x,r,0;t)  -S  ^ J  )         fb(p,q>)- 

k=1  rh  0L 


r[p  +  *  (p)]2  sin(cp+5,-  0) 

•  3— dcp  dp  (9) 

D 

k   iN(nt&k)  fl  p0T 

ut  JC«,r^,t)  -£  "6  4tf  .    J  J    fb(p,cp). 

k=1  rh  ^L 


2   2    A 
(x-  A(p)cp[p  +  A  (p)]2  cos(cp+5,  -  0) 


D' 


-x dcp  dp    (10) 


K    iN(ftt40T-5k)  (l   p00 


k=1  rh  0T 


r[p  +  A  (p)]2  sin(cp+5k  -  0) 

3 dcp  dp  (11) 

D 
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(r) 


u^J(x,r,e,t)  =S 


„=*-* j    J    rr(p).p"^ 

k=1  rh  eT 


(x-  A(p)cp)[p  +  A  (cp)]2    cos(cp+8,   -0) 


dcp   dp  (12) 


t0  *     e1^*-5^    fc**    . 

ua  'C*.«»«,t>  -E  R ^   J      rf(P,cp) 

rh  SL 


[p    -  pr   cos(cp  +  5,  -  0)  ] 


dtp   dp 


(13) 


_        iN(flt-6,)     .Ufi 

rf)  5  e         k    r  r  t 

iO    y(x,r,0,t)    =2J  rr- 

k  =  l 


n 


rf(p,cp) 


C   A(p)[r-  p   cos(cp+5,-  0)] 

.{ . i + 


D 


p[x-A(p)cp]    sin(cp+ 5k-  0)>| 
+  ^ j   dcp   dp 


(HO 


^>(x,r,a,t)  -S  — 

k=1  rh  *T 


J     J      w 


(p)-p 


•iNcp 


[p    -  pr    cos(cp+  5,  -  0)  ] 


dcp   dp 


D" 


K      eiN(^^T-5k)    rl 
u  w  (x,r,e,t)   =E    r— 

k=l 

rh   dT 


n\cp)-p-i% 


(15) 


{ 


A(p)[r-  p    cos(rp+B    -  0)         p[x-A(p)cp]sin(cp+ 5,  -   0) 


+ 


j    d'P    dP 


(16) 
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where 

TK(p,Cp)  =  complex  amplitude  of  the  bound  vorticity 
T   (p)    =  complex  amplitude  of  the  radical  vorticity 

in  the  wake 
Tc(pj(p)  =  complex  amplitude  of  the  streamwise  vorticity 

in  the  blade 

T  (p)    =  complex  amplitude  of  the  streamwise  vorticity 
w 

in  the  wake. 

3.   Velocity  Potential 

As  before  we  separate  the  effects  due  to  source  and  vortex 

distributions. 

(i)  Velocity  Potential  due  to  Source  Distribution 

The  velocity  potential  at  a  point  due  to  a  source  of  unit 

strength  located  at  a  distance,  d,  is  given  by  [16] 

Multiplying  by  the  source  strength,  integrating  over  the  first 
blade  and  summing  over  the  K  blades  gives: 

1 1  K     iN(nt"Sk>  r1  ceT  .«    s 

.W(,,tfe,0  -S  ^Ji j  j     vo(P).p-^  .  |i  • 

k=1  rh  eL 


(17) 
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(ii)  Velocity  Potential  due  to  Vortex  Distribution 
Instead  of  trying  to  find  an  exact  expression  for  the  velocity 
potential,  which  would  be  rather  cumbersome,  we  make  use  of  the  fact 
that  the  unsteady  effects  become  important  only  at  a  distance  from 
the  propeller  and  that  a  reentrant  vortex  is,  in  a  sense,  equivalent 
to  a  uniform  distribution  of  double  sources  over  any  surface  bounded 
by  it.   It  is  shown  in  Appendix  A  that  if  we  substitute  a  reentrant 
vortex  by  an  equivalent  concentrated  doublet,  the  error  in  the  value 
of  the  velocity  potential  decreases  rapidly  with  distance. 

In  this  way  we  can  compute  the  velocity  potential  by  summing 
a  set  of  discrete  values  given  by  expressions  of  the  form  (see  Ap- 
pendix A) 

S  cos  H 


$  = 


kir  D" 


where 


S  =  strength  of  the  concentrated  doublet 
Y  =  angle  between  the  doublet  axis  and  the  line  from 
the  doublet  to  the  field  point. 
The  value  of  S  will  be  found  later  and  cos  ¥  is  easily  computed 
from 

ii-D 


cos  Y  = 


|n|.D 


where 


n  =  vector  direction  of  the  double  axis 
D  =  [x  -  A(p)cp]i+  [r  cosS  -  p  cos(cp  +o\  )  ^  J  + 
+  [r  sin0  -  p  sin(cp+5,)]k 
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Using  (5)  to  obtain  n,  it  follows  that 

[x  -  A(p)cp]cosf3(p)  +  r  sinf3(p)cos(cp  +  5,  -  Q) 
cos  ¥  = (19) 


k .      The  Loading  Modes 

For  the  purposes  of  this  work  it  is  assumed  that  the  distri- 
bution of  the  loading  on  the  propeller  blades  has  already  been  de- 
termined by  any  of  the  known  methods  (examples  [2],  [3],  [7],  [12] 
and  [15])  which  solve  the  unsteady  lifting  surface  problem  of  a 
propeller  operating  in  a  wake.   The  usual  procedure  in  these  methods 
is  to  assume  that  the  loading  may  be  expressed  in  separable  form  as 
the  product  of  a  function  of  the  spanwise  coordinates  and  a  function 
of  the  chordwise  coordinates.   In  selecting  these  functions,  one 
makes  use  of  the  known  results  of  the  theory  for  incompressible  flow 
in  two  dimensions  and  of  the  lifting  line  theory  for  the  higher 
aspect  ratio  wings.   This  leads  to  functions  of  the  type 


Ap  =  Z)  F  (s)  X  G  (r) 

r\        n  n 

n=0 


where 


Ap    =  loading  per  unit  area 

F  (s)  =  chordwise  function 
n 

G  (r)  =  spanwise  function, 
n       r 

Two  dimensional  airfoil  theory  suggests  that  we  choose  F   of 


the  form,  [5] 


i 
Fn(l)  =  [frlT  +  sin[n(cos*|)] 
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where 


I   =  chord  length 
Introducing  a  new  coordinate  9    defined  by 

s   +  sle    i  * 

s  =  r —   "  f  COS0  -  0  -  ¥ 

F  becomes 
n 

F  (0)  =  cotan  -r  +  cos(nTT)  sin(n0) 
n  I 

The  spanwise  function  G  (r)  is  either  determined  directly  with 
the  aid  of  a  specially  selected  interpolation  function  [5]  or  ex- 
pressed in  the  form 

G(r)  =S  A  H  (r) 
n  nm  m 

m 

and  then  determined  by  the  solution  of  the  integral  equation  resulting 
from  the  formulation  of  the  lifting  surface  problem.   In  order  to  ob- 
tain a  numerical  solution  G  (r)  is  usually  approximated  by  a  stepwise 
distribution. 

The  loading  then  will  be  of  the  form 

JJ      00 

Ap(r,9)=  L  (r)  cot  %  +  D  L  (r)  sin(ne)  (20) 

o        z    n   n 
n=l 

where  L   and  L   are  stepwise  functions  of  the  radius.   This  is  called 
on 

the  Birnbaum  distribution. 

When  the  loading  is  expressed  in  terms  of  the  bound  vorticity, 
it  is  necessary  to  add  an  extra  term  to  the  Birnbaum  modes  in  order 
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to  satisfy  the  Kutta  condition  at  the  trailing  edge  [13].   In  this 
case  it  becomes 

JJ      00 

Tb(r,e)  =  AQ(r)  cot  |  +  £  An(r)  sin(n0)  +  f(r,0)      (21) 

n=l 


where  the  Kutta  condition  requires  that 


where 


f(r,Tr)  =  Tb(r,TT)  =  -iknr  (22) 


=  N _  (l  i(r) 


n    2  V  (r) 
o 

r  =  chordwise  integration  of  T, 
According  to  Brown  [3]  a  convenient  choice  for  f(r,0)  is 

r/    £\       -11-     1   -  COS0 

f(r,0)  =  -iknr  x  g 

After  substituting  into  (21)  and  integrating  over  the  chord,  it  fol- 
lows that 

ik 


f(r,e)  =  -7T  x  -  +;k  X  (Ao  +  i  A  )  X  '    ~2C°SU  (23) 


1   -  COS0 
X  <,A   f  f  A  )     X 
O    ^    1 

n 
It  must  be  noted  that  there  is  no  correspondence  (except  for 

the  zeroth  harmonic)  between  the  coefficients  of  (20)  and  (21). 

However,  when  the  loading  is  given  in  terms  of  pressure  distribution 

one  can  always  find  the  corresponding  vorticity  distribution  from 

the  solution  of  the  following  Volterra  equation  of  the  second  kind 

[U] 

Ap(rC°    -£j     rb(r,s,t)ds  +  v0(r)  rbCr,.,t)  (24) 

f  sL 

A. 

in    terms   of  0    the    above   expression   becomes 
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(r^}  =  p   V  (r)  "  ikn(r)  J  Sin^  Tb(r^A)  d^  (25) 


f   o  o 


5.   The  Vorticity  Relations 

The  pulsating  bound  vorticity  is  expressed  in  the  form 

c        iN(fit-5  ) 
Tb(p,cp,k,t)  =  ?b(p,Cp)  e  (26) 

All  the  other  components  of  the  vorticity  can  be  expressed  in 

terms  of  y,,    as  required  to  preserve  continuity  of  vorticity  (see  [1] 

and  [2]  for  detailed  derivation) 

iN(J2t-6k) 
T<:(p/p,k,t)  =  rf(p,cp)  e  (27a) 


f  ^  S   -         2     2     i 

rf(p,cp)  -  -  J  55  (rb(p,cp)-[p  +  A  (P)]2}*p        (27b) 

iN[at-(cp-0  )-5  ] 
Tw(p,cp,k,t)  =  rw(p)  e  K  (28a) 


f„<p)  -  -  ^  (28b) 


iN[fit-(cp-e  )-s  ] 

Tr(p,cp,k,t)  =  rr(p)  e  (29a) 


fr(p)  -  -   2iN^Cp)  ,  (29b) 

[p  +*(p)]2 


where 


rCp)=J    fb(p,cp)-[p2+  x2(P)]*  dp  Oo) 

0L 


-20- 


Notice  that  r  and  r   (streamwise  and  radial  vorticity  in  the  wake) 
w      r  J 

are  expressed  as  the  product  of  a  complex  amplitude  at  the  trailing 
edge  and  a  sinusoidal  time  and  distance  variation  along  the  helical 
surfaces  in  the  streamwise  direction.   This  form  had  already  been 
implied  in  eqs.  (6) ,  (7),  (10)  and  (11). 
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Chapter  III 
NUMERICAL  SOLUTION 

The  following  is  an  outline  of  the  computational  model  used 
to  solve  the  equations  derived  in  the  previous  chapter. 

The  blades  and  wake  are  divided  into  chordwise  panels  and 
these  panels  are  subdivided  by  radial  grid  lines  as  shown  in  Fig.  6 
The  scheme  for  dividing  the  span  is  shown  below. 


d/4   d/2 


FIG.  5  -  SPAN  SUBDIVISION  FOR  NUMERICAL  SOLUTION 


The  interval  is  first  divided  into  MN  equal  spaces.   The  intervals 
at  each  end  are  then  subdivided  into  half  and  quarter  spaces  as 
shown.   The  result  is  MT  =  MN + 4  chordwise  panels.   The  chord  length, 
pitch,  rake,  bound  circulation  strength  and  source  strength  are  as- 
sumed to  be  independent  of  radius  within  each  of  the  panels,  with 
values  matching  those  of  the  continuous  functions  at  the  mid-radius. 
The  continuous  trailing  vortex  sheet  then  becomes  a  set  of  MT+  1 
concentrated  trailing  vortex  lines  at  each  of  the  panel  boundaries. 
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The  bound  circulation  in  the  blade  and  the  radial  circulation 
in  the  wake  are  assumed  concentrated  at  each  of  the  radial  grid 
lines. 

The  above  assumptions  allow  us  to  write  the  integrals  (7)  to 
(12)  as  summations,  as  follows; 

MT  I 
(s),    n    _\     iNflt  ^  ^p   -ilSta(i)    SB(m.i)  r    x/    ,       . .  .  , 
uv  '(x,r,e,t)  =  e     L  L  e   YV    •  — \f —  [x  -  A(m)  ><p(i)  ] 

a  m=l  i=l  L 

K    -iN5,  f^m+l  a 
S  e    M     %  (31) 

MT  I 

ufs)(x,r,0,t)  -  -eiNfitS  S   e"iN:P(i)  .' M^ 
t  m=l  i=l  2 

K   -iN5,  f^1 

S  e    k  •  sin(cp(i)  +  5k  -6)  •  J    p  2fi         (32) 

k-1  p     D 

m 

(b)                                   iNflt            iiTi       GB(m,i) 
u        (x,v,0;t)    -   -e  -r-  Zj    L>       H — 

m=l    1=1 


Pm+1 


K       -iN5,  '  r;mri  d 

S    e  k    •    sin(cp(i)+5,  -  e)    \  -|"  (33) 


Km 


(b),           .      .              iNfit  i?  i       GB(m,i)       ,        ,,    .        .... 
ufc      (x,r,0,t)    =    -e  Jj    Tj       y         •  (x  -  A(m)  'cp(i)) 


MT   I 
m=l    i=l 

Pm+1 


E     e"iN5    .cos(cp(i)  +  5,  -9)    •    \  -§  (34) 

c=l  K  */  J 


P  D 

^m 


Cr^)  iMfl*  MT    °°  iN(0    -cp(i))  ,    . 

u^r)(x,r,0,t)    =    -eMt    -r-  £    S        e  T  •    ^M 

a  m=l    i=I  2 

K        -iNbk  pPm+1   d 

E     e  .sin(cp(i)  +  5,  -  0)    A  -~  (35) 

k=l  k  J  D3 

Mm 
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MT 


(r) ,  n      s  iNfit  tta    v> 

\xK    '(x,r,0,t)    =   -e  E    £     e 

m=l  i=I 


U<V9<i>\  GR(m).  ;( 


x  -  A(m)  -cp(i)) 


K 

S  e 

k=l 


•iNo, 


cos(cp(i)  +  5,  -  0) 


■S 


m+1 

dp 

P     D3 
"m 


(36) 


where 

SB(m,i)  =  amplitude  of  the  equivalent  source  strength, 

concentrated  at  the  i ' th  grid  line  in  the  m' th  panel. 
GB(m,i)  =  amplitude  of  the  equivalent  bound  vorticity,  con- 
centrated at  the  i ' th  grid  line  in  the  m' th  panel. 
GR(m)    =  amplitude  of  the  radial  vorticity  in  the  wake,  con- 
centrated at  the  i ' th  grid  line  in  the  m'th  panel. 
In  the  above  equations  the  source  and  vortex  strengths  are  non- 
dimensionalized  with  respect  to  27TRV  ,  and  the  perturbation  velocities 

are  non-dimensionalized  with  respect  to  V  . 

r        s 

The  integrals  that  appear  in  the  equations  above  are  evaluated 
between  the  limits  of  integration  corresponding  to  the  boundaries 
of  each  panel  and  can  be  determined  analytically  as  follows 


Pm+l 


J" 


P  3 


P    D 

Km 


r  p  cos(cp(i)  +  6,  -  0)  +  (x  -  7\(m)«cp(i))   4-  r' 
{[x-  X(m)q>(i)]  +  r2sin2(cp(i)  +  6R-  0)}-D 


IP 


m+1 


m 


/m+1 


p  -  r  cos(cp(i)  +  5k-  0) 


P    D 


[(x  -  A(m)cp(i))   +  r2sin  (cp(i)  +5,  -0)]-D 


nri-1 


m 


The  contribution  from  the  trailing  vortices  --  eqs.  (13) -(16)- 
raises  the  problem  of  streamwise  integrals  that  cannot  be  evaluated 
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analytically.   It  was  found  by  Kerwin  [18]  that  for  the  zeroth 
harmonic  a  simple  trapezoidal  integration  is  as  efficient  as  more 
sophisticated  quadrature  formulas,  except  for  nearby  elements.   In 
order  to  investigate  the  application  of  this  argument  to  other  har- 
monics, we  carried  out  the  integration  for  a  single  trailing  vortex 
line  located  at  0.75R,  using  a  Gaussian  quadrature  and  a  trapezoidal 
rule.   The  computation  of  the  perturbation  velocities  was  carried 
for  different  points  around  the  propeller  using  a  third  harmonic  of 
the  loading.   We  found  that  the  results  obtained  with  the  two  for- 
mulas differ  in  all  cases  by  less  than  2  percent.   However,  due  to 
the  pulsating  character  of  the  wake  vorticity  it  is  necessary  to  keep 
the  intervals  of  integration  small  compared  with  the  wavelength  of 
the  perturbation.   This  results  in  a  very  large  number  of  computations 
for  each  fortex  line  (about  1100  for  a  field  point  in  the  plane  of 
the  propeller) .   A  better  alternative  was  found  by  substituting  the 
infinite  integrals  by  an  infinite  series  of  finite  integrals  and  then 
applying  the  Euler  transformation  to  speed  up  the  convergency  of  the 
series.   Basicly  the  method  works  as  follows: 
Given  an  infinite  integral  of  the  type 


■I 


\  eiNx  -f(x)  dx 
0 

it  is  transformed  into  the  infinite  series 
(k+l)ir/N 
I=S  eiNx  .f(x)  dx  = 

k=o  Jk  tt/n 

=  2   (-1)k       elNt.f(t-!-k£)  dt 
k=0      Jo  N 
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Now,    applying   the   Euler    Transformation    [21],    the   above    series 
becomes 

oo      (-1)     A     a 
I   =  T,    r— 


*,      a\-E   (-1)-(J) 

m=0  x    ' 


«k+1  o rt  \m/      k-m 


tt/N 
'k  = 


a,    =    \        elNt-f(t  +  k  Tr/N)dt 
J0 


Sample   computations   done   with    this   method    showed   good   convergence 
with   a    fairly  reduced   number    of   terms.      Also,    it   was    found    that   a 
5-point   Gaussian    quadrature  was   adequate    for    the   computation   of   the 
terms   of   the   series. 

When    the   above    transformations   are   applied    to   eqs.     (15)    and 
(16),    we   obtain 

/    v  .MO*.   MT+1       iN0rp(m)      __  f   v      K        -1MB, 

(w).  a      *  iNflt  v  T  TL(m)     <p  k 

u        (x,r,0,t)    =   e  l_j         e  •   — r —    2j     e 

m=1  C         k=1 


oo      (-1)J    .   AJa 

'  ^     1+\ ~   >  (37) 

JO  2J  +  l 


ru/N  -ui(wT) 

a      =    \        e  ■    FWA(£  +  0    +  k  7r/N)d£ 

J        JQ  T 


(w)  /  o    *.\  lNOt       -^  T  TL(m)  -<p  k 

u";   '(x,r,0,t)    =   e  •  Li       e  •   — 5 — 2^     e 

m=1  *     k=1 

00      (-1)J    •   AJa 

*  S     J+1 °  1  <38) 

J=l  2J     ' 


TT/N 

f     -iN(i+e  ) 

a     =    \      e  •    FWT(|+0_+ k  1T/N)dI 

J        JQ  T 
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where 

TL(m)      =   amplitude   of   the   m' th   trailing   vortex,    on   the  wake 

p    [p    -r    cosCcp+S,-©)] 
"mm  k 


FWA(cp)    = 


D3 


FWT(cp)    =    {-p    [A(m)cos(cp+  5,-  0)  -  (x  -  ?\(m) -cp)  sin(cp  + 5,-  0)  ]   + 
r-A(m)}/D3 

It  is  clear  that  eqs.  (37)  and  (38)  cannot  be  applied  to  the  zeroth 
harmonic  of  the  loading  because  the  upper  limit  of  the  integrals 
would  go  to  infinity.   Instead,  we  will  use  the  trapezoidal  rule 
which  was  found  appropriate  for  this  case,  as  mentioned  before. 
Using  a  trapezoidal  rule,  eqs.  (15)  and  (16)  become 

/  \  .„n  MT+1  oo   iN(9  -cp(i))    _  ,       .%,„,/   .  ,N 

(w)^    «  „\     iNfit  v    v       T  ^K  TL(m,i)  +  TL(m,i-l) 

u    (x,r,0,t)  =  e     Li        2j     e  ■  ?    . L-1 - 

m=l   i=I 

K   -iN5,  p  -r  cos(cp(i)+ 5,  -  0) 
p  •  Zj     e      x (39) 

m  k=1  D3 

(w),  a.        mat    M™f  iN(e  -cp(i))  TL(mii)  +  TL(mtl.1) 

u    (x,r,0,t)  =  e     •  Zj  Li     e  •  A — 2 ' 

m=1  i=I  4 

K   ~m\ 

£  e     {-P  [A(m)cos(cp(i)  +  6,  -  0)  - 
k=1  m  R 

(x  -  A(m)cp(i))-sin(cp(i)  +5k  -  0)]   r  A(m)  }/D3 

(40) 
The  contribution  from  the  trailing  vortices  on  the  blade  --  eqs.  (13) 
and  (14)  --  are  computed  using  a  trapezoidal  integration,  for  all  the 
harmonics, 
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MT+1     T 

(f),  .       .  iNfit     iA       fv  GT(m,i)  +  GT(m,i-l) 

ua      (x,T,9,t)    =   e  Li      Li £ ■ 

m=1    i=1 


K        -iNB,    p    -r    cos(cp(i)  +  5,    -  0) 


k  Km  *YX'         k 


•  p       2j     e  - (41) 


m  k=1  D3 


MT-M        T 

(f),           a    „*           iNfit     £      <A     GT(m,i)+GT(m,i-l) 
ufc      (x,r,0,t)    =  e  L      Zj    ' 2 • 

m=1    i=1 

K        ~mbk 

'  Z)     e  C-pm[^(m)cos(cp(i)+S,  -  0)    - 

km  k 

=1 


(x  -  A(m)cp(i)'sin(cp(i)  +  5    -0)]+r   A(m)  }/D3        (42) 


where 


GT(m,i)    =  Amplitude   of   the  m' th   trailing  vortex   located 
on   the  blade ,    between  the   i'th   and   the   i'th+1 
grid   lines. 
The    triple    summations    corresponding   to    the   velocity   potential 
due   to   the    source   and  vortex  distributions    --   eqs.    (17)    and    (18)--   are 

MT         T 

$(s)(x,r,e)t)    -   -e1Nnt     £      S     e-«*P(i).    SBJe.11    . 

m=1    1=1 

k=1  p 

,,%                                    .„r,       MT      I      __,       .v         K      -iN5, 
,(b).                  ,              iNfit     ^      ^     GP(m,i)        ^  k 

$        (x^r^t)    =    -e  Li      Li     — J—2 ?_y  e 

m=1    i=1  k=1 


(x  -  A(m) -cp(i))  cos   (3(m)+r    sin   (3(m)    cos  (cp(i)  +  5,    -  0) 


k 


(44) 
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.*.(b),    „   v      lNflt  <c^   v->  GL(m.i)  v^       k 
0V  y(x,r,9,t)  =  -e      EL  — ^~       S  e 

W  m=1  1-1    *    k=1 


(x  -  A(m) -cp(i))cos  P(m)  +  r  sin  P(m)  cos(<p(i)  +  5  -  9) 


(45) 


where 


<£,    =  velocity  potential  due  to  the  vorticity  in  the  blade. 

3>  =  velocity  potential  due  to  the  vorticity  in  the  wake, 

w 

GP(m,i)  =  amplitude  of  the  concentrate  vortex  equivalent  to 

GBOji)  . 
GL(m;i)  =  amplitude  of  the  concentrate  vortex  equivalent  to 

GR(m,i). 
The  integral  appearing  in  eq.  (43)  is  evaluated  analytically 


Art-1  dp 

J     —  =  to   |2(D+  p-  r  cos(cp(i)  +  5k-  9))| 


rn+1 


m 


The  value  of  the  bound  circulation  is  obtained  either  directly 
if  the  loading  is  expressed  as  in  eq.  (21)  or  after  solving  a  Volterra 
equation  (25)  if  it  is  expressed  as  in  eq.  (20). 

The  solution  of  the  Volterra  equation  can  be  obtained  using  the 
Liouville-Neumann  method  of  successive  substitutions  [20].   This  gives 


Tb(r'e)  =    -V  (r)  ^AP^r^>  "  (iKn(r)  + 


f   o 


£  (-D 

k=1 


k    ^ikn(r)) 


k+1 


k: 


-ik  (r)  r 
k^      n     ' 

cos  9   e        •  \  Ap(r ,  t)  sintdt] 

J0 


(46) 
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When  the  pressure  loading  is  given  in  the  form  of  Birnbaum  modes 
(20),  the  above  integral  is  evaluated  using  the  relations 

\  t  ^       /\ 

\   cotan  -r  sin  t  dt  =  6  -  sin0 

Jo     l 


I 


9  " 

sin  nt  sin  t  dt  =  »  -  f  sin(20)  ,    n  =  1 


sin((n-  1)-0)     sin((n  +1)-9)      , 
2(n-  1)      "    2(n+1)   ~>   n^ 
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Chapter  IV 

RESULTS 

The  numerical  solution  developed  in  the  preceding  chapter  was 
put  into  computer  form;  the  details  of  the  program  are  presented  in 
Appendix  B.   In  order  to  test  the  workability  of  this  program  we 
used  it  to  compute  the  pressure  at  an  upstream  field  point  using  the 
NSRDC  Propeller  No.  14-118  [22].   The  design  characteristics  and  oper- 
ating conditions  of  this  propeller  are  represented  in  Figs.  7  and  8. 
We  considered  both  the  open  water  condition  and  the  nonuniform  wake 
obtained  with  a  screen  which  produces  a  circumf erentially  varying 
sinusoidal  wake  composed  only  of  multiplies  of  blade  harmonics. 
Unfortunately,  we  did  not  have  information  about  the  loading  distri- 
bution for  those  conditions  and  therefore  we  used  an  arbitrary  number 
of  modes  of  the  Birnbaum  distribution  with  coefficients  equal  to 
unit.   In  order  to  test  the  solution  of  the  Volterra  equation  (k6) 
the  input  loading  was  expressed  in  terms  of  pressure  distribution 
rather  than  vorticity  distribution. 

We  made  several  sample  calculations  in  order  to  evaluate  the 
effect  of  changing  the  grid  spacings  and  to  determine  how  far  down- 
stream we  needed  to  carry  the  streamwise  integrations  in  the  wake 
(Table  1  shows  a  sample  of  the  computer  output  for  one  of  the  cal- 
culations).  We  found  that  for  the  field  point  considered  (x  =  -.8, 
r  =  1«2,  Q    =  0)  ;  we  needed  about  20  chordwise  panels  and  had  to  keep 
the  angular  grid  spacing  down  to  2  degrees.   The  contribution  due  to 
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r/R=1.0 


D 

12  Inches 

Dhub~ 

2.4  inche 

s 

r/R 

1/D 

T/l 

P/D 

0.2 

0.3200 

0.1028 

1.0860 

0.3 

0.3635 

0.0776 

1.0845 

0.4 

0.4048 

0.0590 

1.0823 

0.5 

0.4392 

0.0451 

1.0796 

0.6 

0.4610 

0.0348 

1.0780 

0.7 

0.4622 

0.0271 

1.0766 

0.8 

0.4347 

0.0210 

1.0750 

0.9 

0.3613 

0.0166 

1.0735 

1.0 

1.0700 

FIG.  ?  -  PROPELLER  4118  DESIGN  CHARACTERISTICS 
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FIG,  8  -  RADIAL  DISTRIBUTION  OF  THE  AMPLITUDES  AND 
PHASES  OF  THE  THIRD  HARMONIC  OF  THE 
CIRCUNFERENTIAL  WAKE  VARIATION 
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********  PRESSURE  FIELD  AROUND  A  MARINE  PROPELLER  ***-*>.** 


PROPELLER  CHARACTERISTICS 

TYPE         N3RDC,  PROPELLER  NO .  4118 

DIAMETER      12.0000  INCHES 

HUB  DIAMETER    2.1+000  INCHES 

NO.  OF  BLADES   3 

RPM  1080.0000 

BLADE  OUTLINE  (UNITS -FRACTION  OF  RADIUS) 
RADII    LEAD. EDGE  TRAIL. EDGE  MAX. THICK. 


0.250 

-0.3417 

0.3417 

0.0605 

0.300 

-0.3635 

0.3635 

0.0564 

0.400 

-0.4048 

0.4048 

0.0478 

0.500 

-0.^392 

0.4392 

0.0396 

0.600 

-0.4610 

0.4610 

0.0321 

0.700 

-0.4615 

0,4615 

0.0250 

0.800 

-0.4347 

0.4347 

0.0183 

0.900 

-0.3613 

0.3613 

0.0120 

0.950 

-0.2500 

0.2500 

0.0065 

1.000 

0.0 

0.0 

0.0 

GRID  SPACING 

NO.  OP  CHORDWISE  PANELS    20 

ANGULAR  CHORDWISE  SPACING  0.0349066  RADIANS 


TABLE  1  -  SAMPLE  COMPUTER  OUTPUT 
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INFLCW  CHARACTERISTICS 


SYMBOLS  AND  UNITS 

VS  SHIP  SPEED,  KNOTS 

RADII  FRACTION  OF  RADIUS 

XVXO  ZEROTK  HARMONIC  OF  INFLOW  /  VS 

TANBETA  TANGENT  OF  ADVANCE  ANGLE 

NHARM  HARMONIC  NUMBER 

XVXN  AMPLITUDE  OF  NHARM  OF  INFLOW  /  VS 

XPHN  PHASE  OF  NHARM  OF  INFLOW,  RADIANS 


vs= 

8. 8?80 

NHARM=   3 

RADII 

XVXO 

TANBETA 

XVXN 

XPHN 

0.250 

1.0000 

1.0606 

0.0800 

-0.2269 

0.300 

1.0000 

0.8838 

0.1000 

-0.1920 

0.400 

1.0000 

0.6629 

0.1500 

-0.1047 

0.500 

1.0000 

0.5303 

0.2000 

0.0 

0.600 

1.0000 

0.4419 

0.2500 

0.0873 

0.700 

1.0000 

0.3788 

0.2900 

0.0 

0.800 

1.0000 

0.3314 

0.3200 

-0.0314 

0.900 

1.0000 

0.2946 

0.3400 

0.0 

0.950 

1.0000 

0.2791 

0.3450 

0.0436 

TABLE  1  (cont.) 
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LOADING  DISTRIBUTION 
(VORTICITY  MODES) 


RADII 


0.250 
0.250 
0.250 
0.300 
0.300 
0.300 
0.400 
0.400 
0.400 

0.500 
0.500 
0.500 

0.600 

0.600 

0.600 

0.700 
0.700 
0.700 

0.800 
0.800 
0,800 
0.900 
0.900 
0.900 
0.950 
0.950 
0.950 


MODE 


1 
2 

3 

1 
2 

3 

1 
2 

3 

1 
2 

3 
1 
2 

3 
1 
2 

3 
1 
2 

3 
1 
2 

3 
1 
2 
3 


COEF. 


REAL 

0.014860 
0.014860 
0.014860 
0.014860 
0.014860 
0.014860 
0.014860 
0.014860 
0.014860 
0.014860 
0.014860 
0.014860 
0.014860 
0.014860 
0.014860 
0.014860 
0.014860 
0.014860 
0.014860 
0.014860 
0.014860 
0.014860 
0.014860 
0.014860 
0.014860 
0.014860 
0.014860 


IMAG . 

0.014860 
0.014860 
0.014860 
0.014860 
0.014860 
0.014860 
0.014860 
0.014860 
0.014860 
0.014860 
0.014860 
0.014860 
0.014860 
0.014860 
0.014860 
0.014860 
0.014860 
0.014860 
0.014860 
0.014860 
0.014860 
0.014860 
0.014860 
0.014860 
0.014860 
0.014860 
0.014860 


TABLE  1  (CONT) 
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***  PRESSURE  FIELD  RESULTS  *** 


SYMBOLS  AND  UNITS 

EX     AXIAL  COORDINATE  OF  FIELD  POINT, 

FRACTION  OF  RADIUS 
ER     RADIAL  COORDINATE  OF  FIELD  POINT, 

FRACTION  OF  RADIUS 
TETA   ANGULAR  COORDINATE  OF  FIELD  POINT, 

RADIANS 
NT     NONDIMENSIONAL  TIME,  NT=TIME*RPM/60 
KP     PRESSURE  COEFF.,  P/(DENSITY*( (RPM/ 

60)**2)*(DIAMETER**2) ) 


EX=  -0,8000   ER=  1.2000   TETA=  0.0 


1  -  TOTAL  PRESSURE 


NT 


KP 


REAL 


IMAG 


0.0        0.109260  0.064822 

0.0625  -0.002797  0.087227 

0.1250  -0.101459  0.032849 

0.1875  -0.086291  -0,053374 

0.2500     0.115199  -O.I5675I 

0.3125     0.085257  -0.032765 

0.3750     0.021487  0.064308 

0.4375  -0.051879  0.049902 

0.5000  -0.0304070  -0.024493 

O.5625  -0.045830  -0.065644 

0.6250     O.036785  0.004191 

0.6875     0.028781  0.083445 

0.7500  -0.047159  0.067708 

0.8125  -0.08048  -0.061579 

0.8750     0.040268  -0.05525^ 

0.9375     0.117604  0.020541 
TABLE  1  (cont) 
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2 

-  HARMONIC  COMPONENTS 

HARM. NO. 

KP 

REAL 

IMAG 

0 

0.000563 

0.000171 

1 

-0.000115 

0.001645 

2 

0. 008606 

0.000716 

3 

0.069201 

0.047710 

fc 

-0.001695 

-0.002692 

5 

-0.011514 

-0.011500 

6 

-0.010340 

0.001098 

7 

0.000716 

0.001721 

8 

0.010076 

0.001721 

9 

-0.001120 

-0.009159 

10 

-0.000015 

0.000259 

11 

0.000433 

0.000063 

12 

0.019964 

-0.009159 

13 

0.000831 

-0.000075 

14 

-0.000070 

-0.000130 

15 

0.018597 

0.007825 

TABLE   1    (cont.) 
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the  wake  after  about  7  revolutions  was  found  negligible,  which  agrees 
with  calculation  by  Kerwin  [18]  for  the  zeroth  harmonic.   For  field 
points  down  stream  the  wake  integration  is  carried  down  7  revolutions 
beyond  the  axial  position  of  the  field  point. 

It  is  shown  in  Table  2  how  the  number  of  chordwise  panels  and 
the  angular  grid  spacing  affect  the  value  of  the  harmonic  components 
of  the  pressure. 

Figures  9  and  10  show  the  pressure  variation  and  harmonic 
components  at  the  field  point  due  to  the  zeroth  and  third  harmonic 
of  loading  respectively. 

Since  we  used  an  arbitrary  loading  distribution  we  cannot  assess 
the  accuracy  of  the  results.   This  will  have  to  be  done  at  a  later 
date. 
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No.  of  panels 
Grid  spacing 


A 

15 

4  degrees 


B 

15 

2  degrees 


HARM. NO. 

KP 

REAL   * 

IMAG 

0 

0.0 

0.0 

3 

0.081273 

0.054870 

6 

-0.014081 

0.000345 

9 

0.002288   • 

-0.017413 

Kp 

REAL       IMAG 

0.0  0.0 

0.076421  0.053042 

•0.009018  0.000736 

0.001631  -0.012273 


No.  of  panels 
Grid  spacing 


20 

2  degrees 


HARM. NO. 

KP 
REAL   y 

-  IMAG 

0 

0.0 

0.0 

3 

0.069201 

0.047710 

6 

-0.010340 

0.001098 

9 

0.001120   ■ 

"0. 009159 

TABLE  2  -  EFFECT  OF  NUMBER  OF  CHORDWISE  PANELS 
AND  ANGULAR  GRID  SPACING 
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Field  point:   x/R  =0.8 
r/R  ■  1.2 

0  ■  0. 


HARMONIC  COMPONENTS 

HARM. NO. 

KP 

real 

imag. 

0 

-0.059293 

0.0 

1 

-0.000003 

-0.000004 

2 

-0.000002 

-0.000002 

3 

-0.001161 

-0.064603 

4 

-0.000003 

-0.000001 

5 

0.000000 

-0.000001 

6 

-0.000897 

-0.000608 

7 

-0.000006 

0.000004 

8 

0.000000 

0.000000 

9 

-0.000006 

-0.000004 

FIG.  9.b  -  PRESSURE  AT  UPSTREAM  FIELD  POINT  DUE  TO  ZEROTH 
HARMONIC  OF  LOADING  -  HARMONIC  COMPONENTS 
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Fleld  point:   x/R  =0.8 
r/R  =  1.2 
0  =  0. 


.HARMONIC  COMPONENTS 

HARM. NO. 
0 

real 
0.000563 

KP 

imag. 
0.000171 

1 

-0.000115 

0.0016^5 

2 

0.008606 

0.000716 

3 

0.069201 

0.0^7710 

^ 

-0.001695 

-0.002692 

5 

-0.00151^ 

-0.001500 

6 

-0.0103^0 

0.001098 

7 

0.000716 

0.001721 

8 

0,000076 

-0.000898 

9 

-0.001120 

-0.009159 

FIG. 10. b  -  PRESSURE  AT  UPSTREAM  FIELD  POINT  DUE  TO 
THIRD  HARMONIC  OF  LOADING  -  HARMONIC 
COMPONENTS 
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Chapter  V 
CONCLUSIONS  AND  RECOMMENDATIONS 

A  computerized  procedure  has  been  developed  for  the  calculation 
of  the  pressure  field  generated  by  a  propeller  operating  in  an  un- 
steady flow.   We  believe  that  it  represents  an  important  improvement 
over  the  existing  methods  based  on  lifting  line  theory  and  circum- 
ferentially  uniform  wake.   It  is  also  a  useful  complement  to  the 
recently  published  work  by  Jacobs  e_t  aJ.  [19]  which  develops  a 
different  approach  to  the  same  problem. 

A  considerable  amount  of  time  was  put  into  the  optimization  of 
the  program  but  it  is  probable  that  it  can  still  be  ameliorated. 
One  of  the  most  time  consuming  operations  --  the  integration  of  the 
streamwise  vorticity  in  the  wake  --  was  successfully  reduced  by  the 
use  of  the  Euler  Transformation.   The  two  other  operations  whose 
improvement  would  have  an  important  effect  on  the  over -all  improvement 
of  the  program  are  the  integration  of  the  radial  vorticity  in  the 
wake  and  the  calculation  of  the  velocity  potential. 

Sample  calculations  were  made  to  determine  the  sensibility  of 
the  results  to  the  number  of  chordwise  panels,  angular  grid  spacing 
and  how  far  downstream  integrations  were  carried.   Assessment  of  the 
accuracy  of  the  results  was  not  possible  since  we  used  an  arbitrary 
loading  distribution. 

Certainly,  the  next  thing  that  needs  to  be  done  is  to  compare 
results  with  experimental  data  or  with  the  results  obtained  in  [19]. 
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This  requires  that  the  loading  distribution  be  obtained  by  any 
of  the  existing  methods  referred  to  before.   The  program  accepts 
the  input  loading  distribution  in  terms  of  pressure  or  circulation, 
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Appendix  A 
VELOCITY  POTENTIAL  INDUCED  BY  A  SINGLE  VORTEX 

In  this  appendix  we  compare  the  velocity  potential  induced  by 
a  single  reentrant  vortex  with  the  one  induced  by  a  concentrated 
doublet  of  equivalent  strength. 

The  velocity  potential  due  to  a  vortex  is  given  by 


^-iff^*  <*-» 


where 

S  =  vortex  strength 

A  =  area  enclosed  by  the  vortex 

¥  =  angle  between  the  normal  to  A  and  the  line  from  A 

to  the  field  point 
D  =  distance  from  A  to  the  field  point. 
Equation  (A-l)  is  similar  to  the  one  that  gives  the  velocity 
potential  induced  by  a  uniform  distribution  over  A  of  source  doublets 
of  strength  S.   It  is  clear  that  the  strength  of  the  equivalent  con- 
centrated doublet  must  be  equal  to  S  x  A. 

In  order  to  simplify  the  calculations,  we  will  consider  a 
rectangular  vortex  centered  at  the  origin,  as  shown  in  Fig.  A-l. 
The  velocity  potential  at  p  is  found  using  eq.  (A-l) 

a  r\  a 

v      *Vj.    [(x-  I)2  I  (y  -  n)2  +  *¥ 


P(x,y,z) 


FIG.  A-l  -  ORIENTATION  OF  REENTRANT  VORTEX 

The  above  integral  can  be  evaluated  analytically  by  integrating 
directly  in  one  of  the  variables  --  i,    say  --  and  then  making  the 
following  changes  of  variables  before  performing  the  next  integration: 


v  = 


t)  -  y 


7  ?    i 

[(x  ar  +  2  r 


u  =  v  +   (v     +  1) 

2 
t   =  u 

The   final   result   is 


where 


$     =  T"   [(tan"1  A   -    tan^B)    -    (tan1  C    -    tan"1  D)  ] 

v       47r  L  J 


2        2  2 

A   =    (y+a)     +z     -(y-a)[(x+a)     +  (y  -  a)     +z    ] 

z(x  +  a) 


(A-2) 
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B   =    Cy+a)2+z2  -(y  +  a)[(x  +  a)2  +  (y  +  a)2  +  22]i 

z(x  +a) 

C  =  (y-a)  +z  -  (y  -a)[(x  -a)  +(y  -a)  +z  ]~2~ 

z(x  -  a) 

D  m    (y  +  a)2  +  z2  -  (y  +  a)r(x-a)2  +  (y  +  a)2  +  z2r 

z(x  -  a) 

The  velocity  potential  induced  at  p  by  a  concentrated  doublet 

2 
of  strength  4a  S,  located  at  the  origin  is  given  by 

7T(x  +y  +z  )2 

It  is  physically  reasonable  to  predict  that  the  largest  dif- 
ferences between  the  values  given  by  eqs.  (A-l)  and  (A-2)  will  be  found 
at  the  field  points  lying  on  the  z  axis.   Figure  A-2  shows  how  those 
differences  vary  with  the  distance  along  the  z  axis. 

It  is  apparent  from  the  analysis  of  Fig.  A-2  that  the  ratio  of 
the  two  values  approaches  very  rapidly  the  unity.   For  z/a  =  8  (i.e., 
for  a  distance  4  times  greater  than  the  side  of  the  original  vortex) 
the  difference  between  the  two  values  is  less  than  one  percent. 
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Appendix  B 
COMPUTER  PROGRAM 

1.   General 

The  program  presented  in  this  appendix  is  essentially  the 
translation  into  FORTRAN  IV  of  the  numerical  solution  outlined  in 
Chapter  III. 

The  initial  section  of  the  program  in  which  the  blade  and  wake 
are  divided  into  panels  is  taken  from  Professor  Kerwin's  program,  [18] 
Also,  the  scheme  used  for  the  zeroth  harmonic  of  the  loading  is  based 
on  Professor  Kerwin's  work. 

The  program  is  written  for  use  on  the  IBM  System/360  digital 
computer.   For  different  processors  slight  changes  might  have  to  be 
made.   The  amount  of  memory  space  required  is  about  150K  bytes  and 
the  computation  time  for  a  single  field  point,  including  harmonic 
analysis  of  the  pressure  is  about  %  minutes  for  the  zeroth  harmonic 
and  5"  minutes  for  higher  harmonics,  using  a  FORTRAN  G  Compiler. 

The  complex  harmonic  analysis  is  performed  by  the  subroutine 
HARM  described  in  [23].   The  final  output  variables  are  subscripted. 
This  was  done  in  anticipation  of  a  plotting  subroutine  which  can  be 
added  in  the  future. 
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2 .      List  of  Key  Variables 

The  following  is  a  list  of  the  principal  variables  in  the  program, 
arranged  in  alphabetical  order. 

A(M,N)    =  Coefficients  of  the  Birnbaum  modes  of  the  loading  function 
BOUNDA   =  Axial  velocity  of  the  field  point,  induced  by  the  bound 

vortices 
BOUNDT    =  Tangential  velocity  at  the  field  point,  induced  by  the 

bound  vortices 
BWAKEA   =  Axial  velocity  at  the  field  point,  induced  by  the  radial 

vortices  in  the  wake 
BWAKET    -  Tangentical  velocity  at  the  field  point,  induced  by  the 

radial  vortices  in  the  wake 
CHOR(M)   =  Local  chord  at  radius  R(M) 
DHUB     =  Hub  diameter 
DIAM     =  Propeller  diameter 

ER       =  Radial  position  of  the  field  point 
EX       =  Axial  position  of  the  field  point 

G(M)      =  Chordwise  integral  of  the  bound  circulation  at  radius  R(M) 
GB(M,N)   =  Bound  circulation  assumed  concentrated  on  the  N ' th  grid 

line  in  the  M'th  chordwise  panel 
GL(M,N)   =  Wake  circulation  assumed  concentrated  between  the 

KWAKE(M)  +  N  and  KWAKE(M)  +  N  +  1  grid  lines  in  the  M'th  chord- 
wise    panel,  for  computing  PTWK 
GLR(M)    =  Advance  coefficient  at  radius  R(M) 
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GLS       =  Advance  coefficient  based  on  VS 

GP        =   Circulation  in  the  blade,  assumed  concentrated  in  the 

middle  of  each  element  of  the  lattice  for  computing  PTBL 
GR        =   Radial  circulation  in  the  wake 
GT(M,N)    =   Streamwise  circulation  in  the  blade,  between  the  N  and 

N  +  1  grid  lines  in  the  M'th  panel 
IDENT     =   Identification  of  the  propeller 
IGB       =   Parameter  that  describes  the  loading  function.   If  IGB =  0 

the  loading  is  in  terms  of  bound  circulation;  if  1GB  >  0 

it  is  in  terms  of  the  pressure  jump  in  the  blade 
KD(N)      =  Angle  between  the  N'th  blade  and  the  YY  axis 
KN(M)      =  Reduced  frequency  at  radius  R(M) 
KWAKE(M)   =   Identification  of  the  first  grid  line  of  the  wake  ar 

radius  RZ(M) 
MN        =  Number  of  equal  spaces  in  the  division  of  the  span  into 

chordwise  panels 
NBLADE    =  Number  of  blades 
NHARM     =  Harmonic  number 
NL        =  Largest  integer  multiple  of  the  grid  spacing  with  value 

less  or  equal  to  TMIN 
NLE(M)     =   Identification  of  the  first  grid  line  in  the  M'th  panel 
NMAX      =  Number  of  loading  modes 
NR        =  Largest  integer  multiple  of  the  grid  spacing  with  value 

less  or  equal  to  TMAX 

NT        =  Total  number  of  grid  lines  in  the  blade 

o 

NTE(M)     =   Identification  of  the  last  grid  in  the  blade  in  the  M'th  panel 
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NTH       =   Parameter  that  defines  the  spacing  between  grid  lines  as 

NTH  X  2  degrees 
NW(M)     =  Total  number  of  grid  lines  spanned  by  the  blade  in  the 

M'th  panel 
PHN(M)     =   Phase  angle  of  the  NHARM'th  harmonic  of  the  inflow  velocity 
PRS(M)     =   Field  point  pressure  when  the  first  blade  is  at  an  angle 

of  M*T7/(8  NBLADE)  with  fixed  vertical  axis 
PTBL      =  Time  derivative  of  the  velocity  potential  at  the  field 

point,  induced  by  the  circulation  in  the  blade 
PTWK      =  Time  derivative  of  the  velocity  potential  at  the  field 

point,  induced  by  the  circulation  in  the  wake 
PTSC      =  Time  derivative  of  the  velocity  potential  at  the  field 

point,  induced  by  the  source  distribution 
R(M)       =  Midradius  of  the  M'th  chordwise  panel 
RPM       =  Rotation  per  minute 

RZ(M)      =  Radius  at  the  lower  bound  of  the  M'th  panel 
SB(M,N)    =   Source  strength  concentrated  on  the  N'th  grid  line  in 

the  M'th  chordwise  panel 
SL(M)      =  Leading  edge  at  radius  R(M) 
SMA(N)     =   Differential  of  the  thickness  at  the  N'th  grid  line,  as 

percent  of  the  max  thickness  in  each  panel 
SOURCA    =  Axial  velocity  at  the  field  point,  induced  by  the  source 

distribution 
SOURCT    =  Tangential  velocity  at  the  field  point,  induced  by  the 

source  distribution 
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ST(M)      =  Trailing  edge  at  radius  R(M) 

STAR(M)     =  Factor  which  multiplied  by  SMA(N)  gives  the  source  strength 

at  the  N'th  grid  line  in  the  M'th  chordwise  panel 
TETA       =  Angular  coordinate  of  the  field  point 
TK(N)       =  Standard  thickness,  fraction  of  the  max  thickness 
TL(M)      =  Amplitude  of  the  streamwise  vorticity  in  the  wake  at 

radius  RZ(M) 
TLE(M)      =  Angular  coordinate  of  the  leading  edge  at  radius  R(M) 
TMAX       =  Extreme  value  of  TTE (M)  over  all  the  panels 
TMIN       =  Extreme  value  of  TLE (M)  over  all  the  panels 
TRAILA     =  Axial  velocity  at  the  field  point ,  induced  by  the 

streamwise  vorticity  in  the  blade 
TRAILT     =  Tangential  velocity  at  the  field  point,  induced  by  the 

streamwise  vorticity  in  the  blade 
TTE(M)      =  Angular  coordinate  of  the  trailing  edge  at  radius  R(M) 
TWAKEA     =  Axial  velocity  at  the  field  point,  induced  by  the  stream- 
wise  vorticity  in  the  wake 
TWAKET     =  Tangential  velocity  at  the  field  point,  induced  by  the 

streamwise  vorticity  in  the  wake 
VFP        =  Inflow  velocity  at  the  field  point 
VS         =  Ship  speed 

VXO(M)      =  Zeroth  harmonic  of  the  inflow  velocity  at  radius  R(M) 
XGL(N)      =  Advance  coefficient  at  radius  XR(N) 
XKN(N)      =  Reduced  frequency  at  radius  XR(N) 
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XPHN(N)     =  Phase  angle  of  the  NHARM'th  harmonic  of  the  inflow 

velocity  at  radius  XR(N) 
XR(N)       =  Input  radius  (see  Input  Description) 
XSL(N)      =  Leading  edge  at  radius  XR(N) 
XST(N)     =  Trailing  edge  at  radius  XR(N) 
XSTAR(N)    =  Factor  of  proportionality  of  the  source  strength  at 

radius  XR(N) 
XTBETA(N)   =  Tangent  of  the  advance  angle  at  radius  XR(N) 
XTZ(N)      =  Maximum  thickness  at  radius  XR(N) 
XVXN(N)     =  NHARM'th  harmonic  of  the  inflow  velocity  at  radius  XR(N) 
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3.   Input  Description 

The  following  sequence  of  input  cards  is  required  for  the 
operation  of  the  program: 

A   IDENTIFICATION 

B   GENERAL  DESCRIPTION 

C   BLADE  OUTLINE 

D  WAKE  CHARACTERISTICS 

E  LOADING  DISTRIBUTION 


1  card 
1  card 

3  cards 

4  cards 

21  cards  max. 


F  FIELD  POINT  DESCRIPTION    1  card 
A  complete  description  of  the  above  input  cards  is  presented  in 
Table  B-l.   The  input  radii  referred  to  in  that  table  are  schematically 
shown  in  the  figure  below. 


XR(1)    XR(2) 


XR(9)    XR(10) 


\  1 


/ 


rh/R 


« >- 

DX 


1.0 


FIG.    B-l    -   INPUT  RADII    OF  PROPELLER  DATA 


DX   is    the    Input    Spacing,    equal    to   -g-  x   interval    from  hub    to   tip, 
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4.   Program  Listings 
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